#include"2BWTSearch_cpu.h"
#include "OccValue_cpu.h"
#include "mian_util.h"
#include "Match_cpu.h"
#include "integer.h"
#include <thrust/sort.h>

uint64 BWT1Mis(uint3* outSA,
                unsigned int* Queries, unsigned int numQuery,
                BWT* bwt, BWT* revbwt) {
	uint64 numSA = 0;
	for(int queryId =0; queryId < numQuery; queryId++) {
		numSA += cpu_BWT1MismatchSearch( outSA + numSA, Queries, queryId, bwt, revbwt); 
	}
	return numSA;
}

uint64 BWTExact(uint3* outSA,
                unsigned int* Queries, unsigned int numQuery,
                BWT* bwt, BWT* revbwt) {
	uint64 numSA = 0;
	for(int queryId =0; queryId < numQuery; queryId++) {
		numSA += cpu_BWTExactMatch(outSA + numSA, Queries, queryId, bwt, revbwt); 
	}
//	printf("exact numSA:%d\n", numSA);
	return numSA;
	
}

uint64 BWT2MisCaseA(uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {
	uint64 numSA = 0;

	for(int queryId = 0; queryId < numQuery; queryId++) {
		numSA += cpu_BWT2MismatchCaseA(outSA + numSA, Queries, queryId, bwt, revbwt);
	}

	return numSA;
}

uint64 BWT2MisCaseB(uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {
	uint64 numSA = 0;

	for(int queryId = 0; queryId < numQuery; queryId++) {
		numSA += cpu_BWT2MismatchCaseB(outSA + numSA, Queries, queryId, bwt, revbwt);
	}

	return numSA;
}

uint64 BWT2MisCaseC(uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {
	uint64 numSA = 0;

	for(int queryId = 0; queryId < numQuery; queryId++) {
		numSA += cpu_BWT2MismatchCaseC(outSA + numSA, Queries, queryId, bwt, revbwt);
	}

	return numSA;
}

uint64 BWT2MisCaseD(uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {
	uint64 numSA = 0;

	for(int queryId = 0; queryId < numQuery; queryId++) {
		numSA += cpu_BWT2MismatchCaseD(outSA + numSA, Queries, queryId, bwt, revbwt);
	}

	return numSA;
}

uint64 CPUBWT2Mismatch(uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {

	uint64 numSA = 0;
/*	unsigned int occValue1[4] = {0,0,0,0};
	unsigned int occValue2[4] = {0,0,0,0};
	unsigned int count = 0;
	unsigned int i,j;
	printf("2mismatch\n");
	for(i=0; i<2897310526; i++) {
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue,
                               i, revbwt->inverseSa0, occValue1); 
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, i, revbwt->inverseSa0, occValue2);
		for(j=0; j<4; j++) {
			if(occValue1[j] != occValue2[j]){
				printf("index: %u, occValue1: %u, occValue2: %u, character: %d\n", i, occValue1[j], occValue2[j], j);
				count++;
			}
		}
		if(count > 40)
			break;
	}
*/
	numSA += BWTExact(outSA + numSA, Queries, numQuery, bwt, revbwt);
	numSA += BWT1Mis(outSA + numSA, Queries, numQuery, bwt, revbwt);
	numSA += BWT2MisCaseA(outSA + numSA, Queries, numQuery, bwt, revbwt);
	numSA += BWT2MisCaseB(outSA + numSA, Queries, numQuery, bwt, revbwt);
	numSA += BWT2MisCaseC(outSA + numSA, Queries, numQuery, bwt, revbwt);
	numSA += BWT2MisCaseD(outSA + numSA, Queries, numQuery, bwt, revbwt);

	return numSA;
}

uint64 CPU2BWTSearch(uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {
	return CPUBWT2Mismatch(outSA, Queries, numQuery, bwt, revbwt) ;
}



uint64 CPU2BWTSearchTwoStrandLimit(long* numMatch, uint3* outSA,
		unsigned int* Queries, unsigned int* revQueries,
		unsigned int numQuery,
		BWT* bwt, BWT* revbwt)  {

	memset(numMatch, 0, sizeof(long)*numQuery*2);
	//const uint64 numSA1 = CPU2BWTSearch_limit(numMatch, outSA, Queries, numQuery, bwt, revbwt); 
	//const uint64 numSA2 = CPU2BWTSearch_limit(numMatch, outSA + numSA1, revQueries, numQuery, bwt, revbwt);	
	//return numSA1 + numSA2;

	return  CPU2BWTSearch_limit(numMatch, outSA, Queries, 2*numQuery, bwt, revbwt);
}



uint64 CPUBWT2Mismatch_limit(long* numMatch, uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {
	uint64 numSA = 0;
	uint64 numSAQuery = 0;
	uint64 tmp = 0;


	for(unsigned int queryId = 0; queryId < numQuery; queryId++) {
		numSAQuery = 0;
		
		if(numMatch[queryId] != -1) {
			numSAQuery += cpu_BWTExactMatch(outSA + numSA + numSAQuery,
                                        Queries, queryId,
                                        bwt, revbwt, &(numMatch[queryId]));
		} else {
			continue;
		}



		if(numMatch[queryId] != -1) {
			numSAQuery += cpu_BWT1MismatchSearch(outSA + numSA + numSAQuery,
                                        Queries, queryId,
                                        bwt, revbwt, &(numMatch[queryId]));
		} else {
			continue;
		}


		if(numMatch[queryId] != -1) {
			numSAQuery += cpu_BWT2MismatchCaseA(outSA + numSA + numSAQuery,	
					Queries, queryId,  
					bwt, revbwt, &(numMatch[queryId]));
		} else {
			continue;
		}

		if(numMatch[queryId] != -1) {
			numSAQuery += cpu_BWT2MismatchCaseB(outSA + numSA + numSAQuery,	
					Queries, queryId,  
					bwt, revbwt, &(numMatch[queryId]));
		} else {
			continue;
		}


		if(numMatch[queryId] != -1) {
			numSAQuery += cpu_BWT2MismatchCaseC(outSA + numSA + numSAQuery,	
					Queries, queryId,  
					bwt, revbwt, &(numMatch[queryId]));
		} else {
			continue;
		}

		if(numMatch[queryId] != -1) {
			numSAQuery += cpu_BWT2MismatchCaseD(outSA + numSA + numSAQuery,	
					Queries, queryId,  
					bwt, revbwt, &(numMatch[queryId]));
		} else {
			continue;
		}



		//finally, also should be checked after  case D
		if(numMatch[queryId] != -1) {
		} else {
			continue;
		}
		numSA += numSAQuery;
	}

	return numSA;
}


uint64 CPU2BWTSearch_limit(long* numMatch, uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {
	return CPUBWT2Mismatch_limit(numMatch, outSA, Queries, numQuery, bwt, revbwt) ;
}

void CPU2BWTInit(BWT* bwt, BWT* revbwt) {
	assert(sizeof(unsigned long long) == 8);
	generateNumOfBit1Table();
	generateReversionTable();

	initShiftTable();
	initCountValue(bwt);

	printf("CPU2BWTInit done.\n");
}


//return: the number of alignment
uint64 reportCPUResult(uint3* sa, const uint64 numSA) {
	uint64 numAlign = 0;

	for(int i = 0; i < numSA; i++) {
		numAlign += (sa[i].z - sa[i].y + 1);
	}

	printTitle("CPU result");
	printf("#alignment: %lu\n", numAlign);
	printLine();

	return numAlign;
}


/**
  @return the number of queries after preprocessing
 */
uint64 preprocess(uint3* outSA, 
		unsigned int* Queries, unsigned int numQuery, 
		BWT* bwt, BWT* revbwt) {

	PRINT_FUNC_NAME;


	int* queryId = (int*)malloc(sizeof(int)*numQuery);
	uint64* depth = (uint64*)malloc(sizeof(uint64)*numQuery);
	memset(depth, 0, sizeof(uint64)*numQuery);

	for(int i = 0; i < numQuery; i++) {
		resetCount();
		CPU2BWTSearch(outSA, Queries + QUERY_SIZE_IN_WORD*i, 1, 
				bwt, revbwt);
		depth[i] = getCount();
		queryId[i] = i;
	}

	thrust::sort_by_key(depth, depth + numQuery, queryId);



	unsigned int* newQueries = 
		(unsigned int*)malloc(sizeof(unsigned int)*QUERY_SIZE_IN_WORD*numQuery);
	for(int i = 0; i < numQuery; i++) {
		memcpy(newQueries + i*QUERY_SIZE_IN_WORD, 
				Queries + queryId[i]*QUERY_SIZE_IN_WORD, 
				sizeof(unsigned int)*QUERY_SIZE_IN_WORD);
	}
	memcpy(Queries, newQueries, sizeof(unsigned int)*QUERY_SIZE_IN_WORD*(numQuery));


	free(newQueries);
	free(depth);
	printf("preprocessing done.\n");

	return numQuery;
}





